A TRIGONOMETRIC METHOD FOR THE LINEAR STOCHASTIC WAVE 

EQUATION 
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Abstract. A fully discrete approximation of the linear stochastic wave equation driven by additive noise is 

presented. A standard finite element method is used for the spatial discretisation and a stochastic trigonometric 

" " ■ scheme for the temporal approximation. This explicit time integrator allows for en'or bounds independent of the 

^ ^ ' space discretisation and thus do not have a step size restiiction as in the often used Stormer-Verlet-leap-frog scheme. 

^ Moreover it enjoys a trace fonmila as does the exact solution of our problem. These favourable properties are 

^>«-' , demonstrated with numeiical experiments. 
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1. Introduction. We consider the numerical discretisation of the linear stochastic wave 
equation with additive noise 

du-Audt = dW in^x(0,oo), 

M = ina^x(0,oo), (1.1) 

m(-,0) = mo, m(-,0) =vo in^, 



where u = u{x,t), '2) C M'', d= 1,2,3, is a bounded convex domain with polygonal boundary 
CsJ \ dSl, and the dot "•" stands for the time derivative. The stochastic process {W(f)}r>o is an 

^ ■ L2(^)-valued g-Wiener process with respect to a normal filtration {,^t}t>o on a filtered 

probability space (ii, J^,P, {j^,}(>o). The initial data uq and vq are ^o-measurable random 
variables. We will numerically solve this problem with a finite element method in space ifTSl 



fT^ I and a stochastic trigonometric method in time 11] and ['4] (see Section[3]). 

There are many reasons to study stochastic wave equations. Let us mention the motion 
(^ I of a suspended cable under wind loading Q; the motion of a strand of DNA in a liquid 

PsJ ■ IS; or the motion of shock waves on the surface of the sun fSl. All these stochastic partial 

differential equations are of course nonlinear and highly nontrivial. But in order to derive 
efficient numerical schemes, we first look at model problems like ( ILlb . 

The numerical analysis of the stochastic wave equation is only in its beginning in com- 
parison with the numerical analysis of parabolic problems. We refer to |JJ and 11251 for 



j^ , spectral-type (spatial) discretisations of our stochastic partial differential equation and to the 



introduction of flSJ for other types of spatial discretisations. We now comment on works 
dealing with the time discretisation of ( ILll i. Strong convergence estimates for implicit one- 
step methods can be found in flTl , despite the main theme of the paper which is weak con- 
vergence. Both for spatial and temporal approximation the order of convergence is found to 
be somewhat lower than the order of regularity, see Remark l272l below. In 1281 the leap-frog 
scheme is applied to the nonlinear stochastic wave equation with space-time white noise on 
the whole line. A strong convergence rate ff{h^'^) is proved, where h is the step size in both 
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time and space, which is in agreement with the order of regularity in this case. The reason for 
this is that the Green's functions of the continuous and the discrete problems coincide at mesh 
points. A similar trick is also used in ll20l and 11211 to derive an "exact" solver Let us finally 
mention the work lfT4l . where error bounds in the p-th mean for general semilinear stochastic 
evolution equations are presented. The authors consider a Fourier Galerkin discretisation in 
space and the exponential Euler scheme in time. This exponential time integrator (see also 
lfT2l . ifTsl . lfT9l and references therein) is, in the linear case, precisely the one that we use ID. 
The paper is organised as follows. Some preUminaries and the main results from lITSi 
on strong convergence estimates for the finite element approximation of our problem are pre- 
sented in Section |2] The stochastic trigonometric scheme is introduced in Section [3] and a 
convergence analysis is carried out in Section 2] A trace formula for the numerical integra- 
tor is obtained in Section |5] and finally in Section |6] numerical experiments demonstrate the 
efficiency of our discretisation. 

2. A finite element approximation of tlie stocliastic wave equation. Before we can 
state the main result on the finite element approximation of lITSl . we must define the spaces, 
norms and notations we will need. Let U and H be separable Hilbert spaces with norms 
\\-\\u, resp. \\-\\h- ^{U,H) denotes the space of bounded linear operators from U to H and 
^2{U,H) the space of Hilbert-Schmidt operators with norm 

\\TU(u.H)--=(t\\Te£Hy^\ 
k=i 

where {ek}k=i is an orthonormal basis of U. If // = U, then ^{U) = ^{U,U) and HS = 
J^2{U,U). Furthermore, if (n,^,P,{^J,>o) is a filtered probabiUty space, then La{^,H) 
is the space of //-valued square integrable random variables with norm 

\\y\\LAa.H)=nHl]"^- 

Let Q € -^{U) be a self-adjoint, positive semidefinite operator The driving stochastic process 
W{t) in (II. Il l is a //-valued g-Wiener process with respect to the filtration {^t}t>i) and has 
the orthogonal expansion 12 3; , Section 2.1] 

^^(0 = L jfHtyj, (2.1) 

where {(77,e^)}°L[ are eigenpairs of Q with orthonormal eigenvectors and {i3j(f)}7=i are 
real-valued mutually independent standard Brownian motions. It is then possible to define 
the stochastic integral /Q4>(s)dW(s) together with Ito's isometry. 
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2 



""=|'||4>Me'/2|||,(^,,)d5, (2.2) 



where O : [0,°°) — > ^[U ,H) is such that the right side is finite. 

For the stochastic wave equation ( II. Il l, we define U = L2{&) and A = —A with D{A) = 
H^{&) n//(}(^). We assume that the covariance operator QofW satisfies 

||V/3-l)/2gl/2||^^<^ (2.3) 

for some j3 > and with the Hilbert-Schmidt norm defined above. If Q is of trace class, i. e., 
Tr(e) = IIG'/^IIhs <°°, thenj3 = 1. If g = A"', i > 0, then /3 < l+s-d/2. This follows 
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from the asymptotic behaviour of the eigenvalues of A, Xj ^ j^''^ . In particular, if Q = I, then 
j3 < 5 and d = 1. Note that we do not assume that A and Q have a common eigenbasis. 
We will use the spaces H" = D{A"'^) for a e E. The corresponding norm is given by 



V a 



:= \\A"/\\ 



L2(S>) 



L^;^'''p.)i.(. 



7=1 



{S)) 



1/2 



where {(Ay,^j)}°Lj are the eigenpairs of A with orthonormal eigenvectors. We also write 
H" =H"xH"-^ and H = H^^H°xH-K 

We use a standard piecewise linear finite element method for the spatial discretisa- 
tion. Let {5^,} be a quasi-uniform family of triangulations of & with /i/f = diam(/ir), h = 
max/fg^j /ijf, and denote by V/, the space of piecewise linear continuous functions with re- 
spect to ,% which vanish on d!^. Hence, V/, C Hq {&) ~ //' . 

We introduce discrete variants of 11 • 11 « and //": 



IK'/z||/,,a = ||A// v/,||i2(5,), v/, eV/,; H/f = V/, equipped with 



\\h,a^ 



where A/, : V/, -^ V/, is the discrete Laplace operator defined by 

(A/,v/,,w/,)l2(5*) = (Vv/,, Vw/,)ij(g), Vw/, e Yh- 

Denoting the velocity of the solution by M2 := tii := ti, one can rewrite (II. Il l as 

dX(f)=AZ(f)dr + BdW(f), f > 0, 
Z(0)=Xo, 



(2.4) 



where A 



/ 

-A 



,B: 



,^ 



and Xo 



The operator A with D(A) 



H =H X //^ is the generator of a strongly continuous semigroup of bounded linear opera- 
tors E(t) = ef^ on //" = H^ x i/^', in fact, a unitary group. 

Let ^/, : 7/*^ -^ V/, and ^/j : H^ -^ V/, denote the orthogonal projectors onto the finite 
element space V/, C Hq{^) = //', where we recall that V/, is the space of piecewise linear 
continuous functions. The finite element approximation of (ILII ) can then be written as 



or in the abstract form 



dM/,,i(f)+A/,M/,,i(f)df = ^/,dH^(f), t >0, 

"/j,l(0) = M/,,0, M/,,2(0) = V/,,0, 



dXft(f) =AhXhit)dt+^hBdW{t), t > 0, 
X/,(0) =Xh,Q, 



(2.5) 



(2.6) 



where A 



/) 







,Xh 



M/i,2 



andX/,,o: 



M/7,0 

v/,,0 



with M/,,0, v/, e V/,. Again, A/, is the 



generator of a Co-semigroup £/,(?) = e''^'' on V/, x V/,. 

It is known, see, e. g., Q Example 5.8] and I1I8L that under assumption ( 12.31) the linear 
stochastic wave equation ( 12.41 ) has a unique weak solution given by 



X{t)=E{t)Xo+ / £(f-5)BdW(i 



(2.7) 
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with mean-square regularity of order )3, 

Similarly, the unique solution of the finite element problem ( |2.6l l is given by 



Xh{t) = £/,(f)X,,,o + / £,,(r - s)^hBm{s). (2.9) 

We quote the following theorem on the convergence of the spatial approximation. 

Theorem 2. 1 (Theorem 5.1 in HSl). Assume that Q satisfies (12.3b for some j3 G [0,4]. 
LetXo = [mo,vo]^ G hI^ = hI^ x hI^^^, X = [ui,U2]^ andXh = [m/i,i,m/i,2]^ ^^ given by (I2.7l l 
a«c/ ( 12. 9t , respectively. Then the following estimates hold fort > 0, where C{t) is an increas- 
ing function of the time t. 

• Ifuhfi = 3^hUQ, v/,,0 = 3^hVo and fi e [0,3], then 

||«,^i(0-«i(OllL,(a,HO)<C(f)/i^'^{||Xo||i^(^,^P) + ||Ai''^-''ei||Hs}. 

• Ifui,o = MhUQ, v/,,0 = >^/jVo and fi G [1,4], r/ze« 

Remark 2.2. Note that the order of convergence in the position, ij5, is lower than the 
order of regularity, /3, in ( 12.8b . This is a known feature of the finite element method for the 
wave equation, see MSH . The upper limits for j3 are only dictated by the fact that the maximal 
order for piecewise linear approximation is 2; higher regularity will not yield higher rate of 
convergence unless higher order finite elements are used, which can be done of course, see 
l\18]l . Similarly, it is shown in M7\ Theorem 4.1] that the order of convergence of implicit 

one-step temporal approximations is ^(A:™"^'^'^+^' '), where k is the steplength and p is the 



order of the method. Thus, p = I and p ~ 2 for the backward Euler-Maruyama and Crank- 
Nicolson-Maruyama methods, respectively. 

We will also use the following relation between A/, and A, see the proof of Theorem 4.4 
in 1161, 

||A«^,,A-"v|li,(5,) < IMli^(^), «e[-i,i], ve//«=L2(^), (2.10) 

where ^/, is the orthogonal projector ^/, : //*' -^ V/,. 

Finally, we remark that the assumption that & is convex and polygonal guarantees that 
the triangulations can be exactly fitted to ^^ and that we have the elliptic regularity || v||^2(^) < 
ClJAvll^^f^) for V e D{A). This simplifies the error analysis of the finite element method. The 
assumption of quasi-uniformity guarantees that we have an inverse inequality and is only used 
in the proof of the case a G [0, 7] of ( 12. 10b . In particular, it is not needed for the proof of 
Theorem lTTl and not for the case j3 = 1 (trace class noise) in the error analysis in Theorem l4. 1 1 
below. 

3. A stochastic trigonometric method for the discretisation in time. In order to dis- 
cretise efficiently the finite element problem (12.5b , or (12.6b . in time one is often interested in 
using explicit methods with large step sizes. A standard approach for the deterministic case is 
the leap-frog scheme, but unfortunately one has a step-size restriction due to stability issues. 
In the present paper, we will consider a stochastic extension of the trigonometric methods. 
The trigonometric methods are particularly well suited for the numerical discretisation of 
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second-order differential equations with highly oscillatory solutions, see ifTOl Chapter XIII] 
for more details. As stated above, the exact solution of (|2.6l l is found by the variation-of- 
constants formula and given by (|2.9l l. We can write £/,(?) as 



Eh{t) 



Ch{t) 

1/2 



-1/2 



Shit) 



-A,'^Si,{t) Chit) 



(3.1) 



1/2, 



1/2, 



with C/,(r) ~ cos(fA^ ) and 5'/,(f) = sin(fAyj'^). Discretising the stochastic integral in the 
sense of Ito, that is, evaluating the integrand at the left-end point of the interval, leads us 
to the stochastic trigonometric method. We let k be the time step size and U^ = m/, o and 
and obtain the numerical scheme U"+^ = Ehik)U" +Ehik)^hBAW'\ that is. 



f/o 



■ VhO, 






I 
n+l 



\"hhik) 



Chik) 
A]/^Shik) Chik) 






-1/2 



Shik) 



Chik) 



3^hAW'\ 



(3.2) 



where AW" ~ H^(f„+i) — W(f„) denotes the Wiener increments. Here we thus get an approx- 
imation [/" « Uhjitn) of the exact solution of our finite element problem at the discrete times 
t„ = nk. 

Remark 3.1. The stochastic trigonometric methods ( 13. 2t are easily adapted to the 
numerical time discretisation of (N -dimensional) systems of nonlinear stochastic differential 
equations of the form 

X{t) + a'xit) = GiXit)) + Wit), 

where (0 £ M.^^^ is a symmetric positive definite matrix and G(x) £ M^ is a smooth nonlin- 
earity. In this case, one obtains the following explicit numerical scheme ^ 



x'i+\ 


_ 


1- 



cos(^a)) (O 'sin(A;a)) 
a)sin(^a)) co^ikco) 



XI 



|('PoG(^^r) +^iG(4>Xf+i)) 



+ 



CO sin(A:a)) 
cos(A:a)) 



(3.3) 



AW", 



where k denotes the step size and AW" = W (?„+ 1 ) — W (r„ ) the Wiener increments. Here ^' = 
Xffikco) and <I> = (pikca), where the filter functions 1/^,0 are even, real-valued functions with 
V/(0) = 0(0) = 1. Moreover, we have ^o = Woikco), ^i = V^i ikco) with even functions i/Aq, i/Ai 
satisfying 1/0(0) — V^i(O) = 1. The purpose of these filter functions is to attenuate numerical 
resonances. Moreover, the choice of the filter functions may also have a substantial infiuence 
on the long-time properties of the method, see l[10\ Chapter XIII] for the deterministic case. 
We will not deal with these issues in the present paper 

Numerical experiments for the nonlinear stochastic wave equation 

du - AMdf = Giu) At + dW 

with a smooth nonlinearity G will be provided in Section |6] in order to demonstrate the effi- 
ciency of this approach. We leave a theoretical investigation of the nonlinear case for future 
works. 

For a more detailed derivation of the trigonometric method and its use for nonlinear wave 
equations we refer to ifTOl Chapter XIII] and ID for the deterministic case and to 121 and H 
for the stochastic case. 
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In the next section we will see that this explicit numerical method permits the use of large 
time step sizes k and that the error bounds are independent of the spatial mesh size h; some 
of these properties are not shared by, for example, the backward Euler-Maruyama scheme, 
the Stormer-Verlet scheme or the Crank-Nicolson-Maruyama scheme, as we will see in the 
numerical experiments in Section |6] 

4. Mean-square convergence analysis. In this section, we will derive mean-square 
error bounds for the stochastic trigonometric method (13.2) . Our main result is a global error 
estimate for the time discretisation in Theorem l4.1l Its proof is based on bounds for the local 
errors in Lemma l4~2l Finally, we formulate an error estimate for the full discretisation. 

Theorem 4.1. Consider the numerical discretisation of (12. 5t by the stochastic trigono- 
metric scheme (|3.2| l with temporal step size k. The global strong errors of the numerical 
scheme satisfy the following estimates: 

• //||A''^"''''2e'''^l|HS <°° for some ^ > 0, then 

• If\\A.^f^-^^I^Q^/-\\m < °° for some jS > 1, then 

The constant C = C{T) is independent of h, k, and n with t„ ^ nk < T. 

For the proof of the above theorem, we will need the following lemma; 



Lemma 4.2. Let the local defects d" = [d",d'2]'^ be defined by 



d'l := T"^' A, "hh{tn+i - 5)^/,dW(5) - A, "hh{k):^hm'\ 

d'l := /'"^' C/,(f„+i - s).'^hAW{s) - Ch{k)^,AW'\ 
Jt„ 

We have the following estimates: 

• //||A('^-l)/2gl/2||^^ <oo for some ^ > 0, then 

E[||<||2^(^)]+E[||A;'/24'||i^(^)]<C^--{2/5+l,3}||^(/3-l)/2gl/2||2^^^ 

• //||A(/^-l)/2ei/2||jj5 < oo for some j3 > 1, then 

E[||A;/V/||i^(^)]+E[||4||2^(^,]<C^--{2/3-1.3}||^(/3-l)/2gl/2|j2^^^ 

The constant C ~ C{T) is independent of h, k, and n with t„ = nk < T. 

Proof. We begin by showing, recall that H^^ = V/, with norm j| -jl/i.o = || ■||l2(®)' 

||(5,(0-5,(.))A,^'''|l^(^0)<C|f-.|^ j3e[0,l]. (4.1) 

For /3 = and v/, e V/, we use the triangle inequaUty and the boundedness of 5'/,(f ): 

||(5/,(f)-5/,(5))v/,||i,(5,) <2|jv/,|ii2(5?) =2||v/,||/,,o. 
For j8 = 1 and v/, e V/, we use the fact that 

iSh{t) - Sh(s))vh = j li,-Sh{r)vhdr ^ j Q(r)Ay^,,dr 
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and hence 

1 /2 

||(5/,(f)-5/,(.?))v/,||i2(s?) < \t-s\\\A,; Vh\\L2m = k-'5|lk/'IU,i- 
A well-known interpolation argument, see e.g. the proof of Theorem 3.5 in IZTl . then yields 

||(5/,(f)-5/,(i)K||i^(^)<c|f-i|'^||v,,||,,^, y,,ey,„ j3e[o,i], 

which is (14. 11 1. 

We now consider t/" with j3 G [0, 1]. By Ito's isometry (I2.2l l and (14. Il l we have 



E[||rf 



"1^2(5*)] 



E 



'n+l 



.-!/■ 



A,, ''-(5,,(Wi -.9)-5/,(A:))^/,dW(.) 



||A,;^/'(5„(.)-5„(fc))^„e'/2||2^^d. 



Liin 



Using also (ITTOl i with a = (j3 - l)/2 e [-^,0] we obtain 

^/A-(^-')/^ll^(H0)iiA(^-')/2ei/2|| 



HS 



< IIA 



(/3-l)/2 



HS 



<c||a(/^-'VV"| 



HS- 



This proves 



Em\\lm]<Ck^P+'\\A^P-'y'Q'/Ym, 



which is the desired bound when j3 G [0, 1]. When/3 > 1, we simply observe that ||A "^ 'II iff//") ^ 
C, so that by the already proven case 



E\\\d' 



<C [\s-kfds\\3^,Q'/YHS<CkV/YHS 
Jo 

<Cfc^^||A(^-l)/2el/2||^^||A-(/^-l)/2||2 



HS 



i^{H") 



<C;t3||A(/3-l)/2gl/2|j2 



HS- 



This is the desired result for j3 > 1 . 

Similarly we find for the second component lij with j3 £ [1,2]: 

E[||41lL(^)]</ll(QW-C,,(fc))A-'^-')/2|l^(^0)d.||Ai'^-''/2^,eV2|,2^^^ 
where, similar to (14.1b . 
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Hence, using also ( 12.101 ) now with a = (j3 — l)/2 e [0, 4], we obtain 



d'i\\U,]<Ck^^-'\\A^^-''l'^Qni 



HS 



for j3 e [1 , 2] . For j3 > 2 the defect is of the order /t^ . 

The bounds for E[||A^j d'lW't f®)] ^^^ ■'^tll^/i '^iW'l f@)] ^^ proved in the same way. D 

We now turn to the proof of our main result on the strong convergence of the numerical 
method ( |33] |. 

Proof. [Proof of Theorem gl] We define F'^ := U'j - Uhj{t„), j = 1,2, and F" 

1 ''^2 J 



[F,", R']^. First of all we remark that 



||f/r-M/,,i(f«)lii2(a//0) = \\^"\\l2{a.H«) -^[\\^"\\l2{9) 



Substituting the exact solution X/, = [m/, i,m/, 2]^ of (I2.6l l into the numerical scheme ( I3.2l i. we 
obtain 



with the defects d" := [d",d2]^ defined in Lemma W2\ and £/,(?) defined in ( 13. It . We thus 
obtain the following formula for the error F"^^: 



F"+'^Ehik)F"+d"=Ehit„+i)F''+Y^Eh{t„^j)dJ = Y,Eh{tn-j)dJ 

7=0 ;=o 



since F" = 0. Taking expectations gives us for the first component 



E IIF, 



riiLjis-)-! 



:E 



:E 



E (C/,(f„-i-,)^{ + A, '^^Shit„_i_j)d^ 



n-l 



Liim 



n-l . "-1 \ 



1-1 



n-l 



"Y^ A-''hh{t„^X-j)d2'Y Ch{tn-l-^)d[] 
j=0 i=0 / 

E \'^^Shit„_i_j)di"Y A/;'/'5/,(r„_i_,)4 

)=0 i=0 
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Here we use the independence of d[ 2 and dl ^ with / , y = 0, . . . , « — 1 for / / 7 to get 



n-\ 



^\\\nfL,(oj^\ =E E(C,,(r„-i-,,K,C,(f„_i_;)«?/) 



n-l 



j=0 
n-l 

+ L (A/;'/'5,,(f„-i-,)a?i',Q(f„-i-;)fif/) 

n-l 
n-l 

<2 



Now we can apply Lemma l4~2l for the estimates of the defects dj and lij and get 



<C(y^)^min{2p,2}|J^(/5-l)/2gl/2|j 



2 
HS- 



Therefore we obtain 



\\U{'-i^„Atn)\Lia.HO) = v/E[||Ff||2^,^)] <C^--{^-'}||A(^-i)/2gi/2||^^ 



for /3 > 0. 

For the second component of F" we obtain 



n-l 



KWFiWli^)] = E ^ (-A,f 5,,(f„-i-,)«''/ + C/,(f„-i-;)«?^' 



7=0 



2 

L2{S*) 



Fl-1 



jiii2 






Thus we get with Lemma |4~2l if j3 > 1 : 



J=o 



<C^mm{2/3-2,2}||^(/3-l)/2g 



1/2||2 



HS 



and 



l|t/2"-«/,,2(r«)IL,(a//0) = ^nmWlirjH < C^™"{'^->-'}||A(P-')/2e>/2| 



HS- 
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D 

We can now collect the convergence results for the space discretisation and for the time 
discretisation. This gives us the following theorem. 

Theorem 4.3. Consider the numerical solution of dl.ll l by the finite element method in 
space with a maximal mesh size h and the numerical scheme ( 13. 2t with a time step size k on 
the time interval [0, T]. Let us denote the discrete time by tn = nk. Let Xq = [uojVq]^ and let 
X = [u\,U2\^ andXh = [m/i.i,m/i.2]^ be siven by (I2.7l i and ( |2.9l l, respectively. t/W^QWi^iQ^uP) ^ 
oo, the following estimates hold for t > 0, where C{t) is an increasing function of the time t. 

• Ifuh,o = ^hUQ, v/,,0 = B^hVQ and ;/ 1| A^'^^^'/^G'/^Hhs < °°for some j3 € [0, 3], then 

• IfuhQ =^/,Mo, V/,,0 = ^hVo and if \\A''-^^^'i I^Q^ I^Was < ^ for some j3 G [1,4], then 

Proof. This follows from Theorems 12. 1 l and 14. 1 I bv the triangle inequality. D 



5. A trace formula for the numerical solution. In this section, we look at a geometric 
property of the exact solution of the wave equation. It is known that, in the deterministic 
setting, the linear wave equation is a Hamiltonian partial differential equation, wherein the 
total energy (or Hamiltonian) of the problem is conserved for all times. However, in the 
stochastic case considered here, the expected value of the energy grows linearly with the time 
t. This is stated in the next theorem for the semidiscretisation of our linear stochastic wave 
equation dl.ll ). For a nonlinear version of this so-called trace formula we refer to [|25l . 

Theorem 5.1. Consider the numerical solution of jl.lb by the finite element method 
in space with a maximal mesh size h. Let X/, = [m/i,i,m/i,2]^ be given by (12.9b . The expected 
value of the energy of the exact solution of the semidiscrete problem ( 12.5b with initial values 
X/,(0) ^ [uh.o.VhflY e L2{^,Vh) satisfies: 



E 



-(llA,yVi(0lli2(5*) + ll«M(0lli2(5*))] =E[^(i|A/yVoiii2(f?) + lK'/.oiiL(5»)) 



l.(||A'/2....|2 ,,,.,..,2 

ifTr(^,e^,) 



for all times t >0. 

Proof. We recall that the solution of ( 12.5b . ^/i(f) = [uh,i{t),uij^2{t)V ^ with initial values 
Xfi{Q) = [uh.o,Vh,o]^ can be written as 



X,it)=Ei,{t)X,,iO)+ f'E,,{t-s)^,,BdW{s). 
Jo 
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Therefore we get for the first summand of the energy, i. e., the potential energy, 
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E 



,1/2 



V «/-i(0IIl2(5*) 



= E 
= E 



1/2, 



2 



A;/"C,,(f)MA.o + 5,,(f )v,,o + j Sh{t - s)^hm{s) 

1/2 9 9 

|A/, Ch(t)uhx)\\l^_^g) + ||5/,(f)v/,_o||i,(<2j) 
/ Sh{t-s),^hAW{s) +2(AyX/,(r)M/,,o,5/,(r)v/,,o) 

Jo L2{9) 

+ 2(^A]/^Ch{tH,oJ'Sh{t-s)^,,dW{s)'^ 
+ 2(s,,{t)vh,o, ['Shit~s)^i,dWis) 



:E 



1/2, 



■h C/,(f)M/,,o|ii,(5,) + |i5/,(f)v/j,o|li2(<J,) 

f Sh{t-s)^hAW{s) \ +2(AyV,,(OM/,,o,5,,(f)v/,,o) 

Jo i2(5*) 



using the fact that the above Ito integrals are normally distributed with mean 0. 
For the second summand we obtain 



E 



|M/.,2(f)llL2(®) 



1/2 

V, 5/,(f)M/,^0llL2(®) + \\Ch{t)Vh.0\\L2{&) 

2 .1/2 



Ci,{t-s),ndWis) 



Lim 



2 (A,/ Ch{t)uh^Q,Sh{t)vhfi) 



Now, we use Ito's isometry to compute, for example. 



E 



Jo L2(3>)i Jo 



Then, combining these expressions and using a trigonometric identity leads to the statement 
of the theorem: 



E 



1 



1/2 



A,VXl(0llt(5*) + ll«/,2(0ilt(S*)) 



:E 



(||A/, l^hfi\\i^(s,) + \\uhfi\\l^(g)) 

^\t\\^hQ"Yns 

-(||A//-M/,^0|lL(5S) + |lM/,.0|li2(i^)) 



+ -fTr(^/,e^/,). 

The last equality follows from the definitions of the HS-norm, of the operator Q and of the 
projector ^/,; 

ll^/,e'/'ll?is =Tr((^;,ei/2)(^;,ei/2)*) = Tr(^;,e^;,). 

This concludes the proof. D 

Remark 5.2. We would like to point out, that an alternative proof of the above result 
can be obtained using ltd' s formula, see for example ^ Theorem 4.17], to the function 



1 



^(f//0 = ^(iiA;/-c/.,iiii(^)+iic/,.2iii,(^)). 
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We are now able to show that the numerical solution given by our stochastic trigonometric 
scheme preserves this geometric property of the exact solution of the finite element problem 
(lOT l. 

Theorem 5.3. Under the assumptions ofTheorem \5.1\ the numerical solution of ( 12.51 1 
by the stochastic trigonometric method ( 13.21 ) with a step size k preserves the linear drift of the 
expected value of the energy, i. e.. 



E 



^(IIA 



h ^riiizi®)' 



^2|Il2(5*)) 



1 



2(11^7 «"^0|lL(5.) + llv'/,,0|li,(^)) 

7f„Tr(^/,e^„) 



for all times t„ ~ nk > 0. 

Proof. The stochastic part of the method can be written as an Ito integral and we obtain 
due to the Ito isometry 



E 



\s,,{k)i^,AW 



n-l||2 



= E 



S„{k)^i,dW{s) 



'n-l 



Lim 



'n-l 



\Sh{k)^hQ 



'/^"^sd.. 



Similarly to the proof of Theorem l5.1l we thus get 



E 



1/2, 



■{w,ru" 



U\L2(S))^ 



|f/^'H2 



L2(5*)^ 



:E 



rl 



l/2,rn-l||2 



{\\^\ru 



\M@)^ 



\u. 



n-l||2 



lL2(5*)^ 



-Tr(^,e^,). 



A recursion now concludes the proof. D 

To conclude this section, we would like to remark that already for stochastic ordinary 
differential equations, the growth rate of the expected energy along the numerical solutions 
given by the forward (or backward) Euler-Maruyama scheme and the midpoint rule, see JS] 
and references therein, is not correct. Indeed, for the forward Euler-Maruyama scheme, one 
has an exponential drift in the expected value of the energy. 

6. Numerical examples. Let us consider the example given in ifTSl : 



dti-AMdr = dW, 

M(0,r) = M(l,r)=0, 

m(x,0) = cos(7r(.\:— 1/2)), m(x,0) = 0, 



(x,r)e (0,1) X (0,1), 

fe(0,l), 

xe(0,l). 



(6.1) 



The solution of this stochastic partial differential equation will now be numerically approx- 
imated with a finite element method in space and the stochastic trigonometric method (13. 2t 
in time. For the below numerical experiments, we will consider two kinds of noise: a space- 
time white noise with covariance operator Q = I and a correlated one. For correlated noise 
we choose Q = A^* with s eM. and recall the relation j5 < I +s — d/2, where d ^ I is the 
dimension of the problem, see the discussion after (12.3) . 

Before we start with our numerical experiments, let us briefly explain how we approxi- 
mate the noise present in the above stochastic partial differential equation. From the Fourier 
expansion ( 12. It , we have for all ;t G V/,: 



i^HAw",x)L,i9) = L Yr^nn^j'XUis 



i&), 
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k=2 , M=1 00 samples 



ior= 



10"' 




-p<1/2 

- Order 1/3 
-p<1 

- OrrJer 2/3 

- p < 3/2 

- Orrier 1 

- p<3 

- Orrier 2 



Figure 6.1. Spatial errors: The Li-error in the first component decreases with order h"!!^. 

where {7;,ey}~Lj are the eigenpairs of the covariance operator Q with orthonormal eigenvec- 
tors {ey}"Lj, and {/3j}°Li are mutually independent standard real-valued Brownian motions 

with Gaussian increments Aj3" = /3;(f„) — j3j(f„_i) ~ y/kJ^{Q, 1). As explained in ifTSl . un- 
der some assumptions on the triangulation and the operator Q, one can approximate the above 
expansion with 



with an integer J >Ni,, where A'/, = dim(V/,), while retaining the convergence rate, to obtain 
the semidiscrete solution, see i 



Xlit) = Ei,{t)Xh,o 



'£rfJ^E,,{t-s)^,Bejdliji 



Figure 16.11 confirms the results on the spatial discretisation of our linear stochastic wave 
equation stated in Theorem 12. 11 The spatial errors in the first component of our problem 
are displayed for various values of the parameter s. On the one hand we consider a space- 
time white noise with Q ~ I, and hence j3 < 1/2, and on the other hand, different corre- 
lated noises with Q ~ A^\ i.e., j3 < l/2 + s. The corresponding convergence rates are 
observed. Here, we simulate the exact solution with the numerical one using a very small 
step size, i.e., ^exact = hexact — 2^^. The expected values are approximated by computing 
averages over M — 100 samples. All the numerical experiments were performed in Matlab 
using specially designed software and the random numbers were generated with the command 
randnC state' ,100). 

We are now interested in the time-discretisation of the above stochastic wave equation 
for various spatial meshes. Figure |62] displays the strong error at time f = 1 in the first 
component of the solution for space-time white noise with s = and for correlated noise with 
i = 1/2, respectively. One observes the order of convergence stated in Theorem 14. H and the 
fact that these errors are independent of the spatial discretisation. Again, the exact solution is 
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s=0, M=1 00 samples 




Figure 6.2. Temporal errors: The Li-error in the first component decreases with order k" and is independent 
of the mesh-grid h. 



approximated by the stochastic trigonometric method with a very small step size ^exact = 2^^. 
We use /zexact — 2^^, 2^'", resp., 2^'^ for the spatial discretisations. Again M ~ 100 samples 
are used for the approximation of the expected values. 

Next, we compare our time integrator with the following classical numerical schemes for 
stochastic differential equations. When applied to the wave equation in the form (12.4b . these 
schemes are: 

1 . The backward Euler-Maruyama scheme X"+ ' = X" + kAX"^ ' + BAW", see for ex- 
ample ifTSl or II22I . The strong rate of convergence for this method is i^(^'"™("'^'''), 
see IITtI Theorem 4.12]. 
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2. A stochastic version of the Stormer-Verlet scheme, writing X — [Xi ,^2]^, 



-1/2 



k 



X^^^'^ =X^ + -AX," +H'(f„+i/2) - W{t„), 

Xf+i =X'^ + kX^'^^^^, 

X^'=xf'/V^AX,"+i+W(Wi)-W(r„+i/2). 

For an application of this scheme to the Langevin equation, we refer to |I24| . We 
were not able to find any references on the strong rate of convergence of this numer- 
ical method. 
3. The Crank-Nicolson-Maruyama scheme ifTTI 

X"+' =X" + '^A{X"+^+X")+BAW". 

The strong rate of convergence is i^(A:™"P"/^'''), see ifTTl Theorem 4. 12]. 
We apply these schemes to the finite element approximation of the linear problem (16. H 
with truncated noise. Note that both the backward Euler-Maruyama scheme and the Crank- 
Nicolson-Maruyama scheme are implicit. Figure l63] presents the various strong convergence 
rates of the above numerical integrators, once with white noise and once with correlated 
noise with Q ~ A^*'^. One observes that the numerical solution given by the Stormer-Verlet 
method explodes for larger values of the step-size k (this computation was stopped when 
the deterministic non-stable regime of the scheme was attained). For all the experiments we 
use /jexact = 2^"' for the spatial discretisation. The reference solution is computed using the 
stochastic trigonometric method with the step size A:exact = 2^'^. Again M = 100 samples are 
used. 

In the following numerical experiment, we are concerned with the trace formula of Sec- 
tion |5] Figure 16.41 illustrates the trace formula of the numerical solution. Here, we choose 
5 = 1/2 and hence )3 < 1 and display the expected value of the energy along the numerical so- 
lution of the above stochastic linear wave equation with mesh grids li — Q.l and ^ = 0. 1 on the 
long time interval [0,500]. We tookM = 15000 samples to approximate the expected energy 
of our problem. A comparison with other time integrators is presented in Figure 16.51 One 
notes that all these numerical schemes do not reproduce the linear growth of the expected en- 
ergy correctly. This fact is already known for the backward Euler-Maruyama scheme applied 
to a finite-dimensional linear stochastic oscillator |!26l. 

Finally we consider a nonlinear stochastic wave equation, the Sine-Gordon equation 
driven by additive noise: 

du-Audt = -sm{u)dt + dW, {x,t) e (0, l)x(0, 1), 

m(x,0)=0, m(-*:,0) = 1j, 3]W, xe(0,l), 

where l/(x) denotes the indicator function for the interval /. The corresponding deterministic 
problem is studied for example in 18j. We solve this problem again with a finite element 
method in space and in time we use the stochastic trigonometric method (13.31) with G{X{t)) = 
— sm{X{t)) and the filter functions proposed in Q: 

V/(^) = sinc3(^), (/)(^)=sinc(^), v/o(^) = cos(^)sinc2(^), xif, {^) = smc^{^), 

where sinc((^) = sm{£,)/£,. In the upper plot of Figure 16.61 we show the expected energy of 
the numerical solution of the Sine-Gordon equation where the covariance operator is given by 
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- Error SV 
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-Order 1/3 
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s=1/2, h=2-^", M=100 samples 



10^' 




- Error BEM 

- Error SV 

- Error CNM 

- Error STM 
-Order 1/2 

- Order 2/3 

- Order 1 



10 10 
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Figure 6.3. L2-error in the first component of the numerical solutions given by the Stormer-Verlet method 
(SV), the backward Euler-Maruyama scheme (BEM), the Crank-Nicolson-Maruyama scheme (CNM) and the 
stochastic trigonometric method (STM). 



Q = I. Even for a large step-size A: = 0. 1 , one can observe the good behaviour of the numerical 
scheme. In the lower figure, we display the convergence rate for the first component with a 
CO variance operator Q = A^^ Again, we approximate the exact solution with a finite element 
solution and the stochastic trigonometric scheme using A;exact = 2^^ and /zexact = 2^^. 
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